Functions and Data Load

loadPkg = function(pkg, character.only = FALSE) { 
  if (!character.only) { pkg <- as.character(substitute(pkg)) }
  if (!require(pkg,character.only=T, quietly =T)) {  install.packages(pkg,dep=T,repos="http://cran.us.r-project.org"); if(!require(pkg,character.only=T)) stop("Package not found") } 
}
loadPkg(knitr)
# Fix outliers

outlierKD2 <- function(df, var, rm=FALSE) { 
    #' Original outlierKD functino by By Klodian Dhana,
    #' https://www.r-bloggers.com/identify-describe-plot-and-remove-the-outliers-from-the-dataset/
    #' Modified to have third argument for removing outliers inwtead of interactive prompt, 
    #' and after removing outlier, original df will not be changed. The function returns the new df, 
    #' which can be saved as original df name if desired.
    #' Check outliers, and option to remove them, save as a new dataframe. 
    #' @param df The dataframe.
    #' @param var The variable in the dataframe to be checked for outliers
    #' @param rm Boolean. Whether to remove outliers or not.
    #' @return The dataframe with outliers replaced by NA if rm==TRUE, or df if nothing changed
    #' @examples
    #' outlierKD2(mydf, height, FALSE)
    #' mydf = outlierKD2(mydf, height, TRUE)
    #' mydfnew = outlierKD2(mydf, height, TRUE)
    dt = df # duplicate the dataframe for potential alteration
    var_name <- eval(substitute(var),eval(dt))
    na1 <- sum(is.na(var_name))
    m1 <- mean(var_name, na.rm = T)
    par(mfrow=c(2, 2), oma=c(0,0,3,0))
    boxplot(var_name, main="With outliers")
    hist(var_name, main="With outliers", xlab=NA, ylab=NA)
    outlier <- boxplot.stats(var_name)$out
    mo <- mean(outlier)
    var_name <- ifelse(var_name %in% outlier, NA, var_name)
    boxplot(var_name, main="Without outliers")
    hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
    title("Outlier Check", outer=TRUE)
    na2 <- sum(is.na(var_name))
    cat("Outliers identified:", na2 - na1, "\n")
    cat("Propotion (%) of outliers:", round((na2 - na1) / sum(!is.na(var_name))*100, 1), "\n")
    cat("Mean of the outliers:", round(mo, 2), "\n")
    m2 <- mean(var_name, na.rm = T)
    cat("Mean without removing outliers:", round(m1, 2), "\n")
    cat("Mean if we remove outliers:", round(m2, 2), "\n")
    
    # response <- readline(prompt="Do you want to remove outliers and to replace with NA? [yes/no]: ")
    # if(response == "y" | response == "yes"){
    if(rm){
        dt[as.character(substitute(var))] <- invisible(var_name)
        #assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
        cat("Outliers successfully removed", "\n")
        return(invisible(dt))
    } else {
        cat("Nothing changed", "\n")
        return(invisible(df))
    }
}
loadPkg('dplyr')
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
loadPkg('readxl')
loadPkg('readr')

#Read in general census data
census_2015 <- read.csv(file = 'acs2015updated.csv')

#Read in Specifically Education Data
edu_data_2015 <- read_csv("ACS_15_5YR_B15003_renamed.csv", skip = 1) %>% mutate(Id2 = as.numeric(Id2))
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   Id = col_character(),
##   Geography = col_character()
## )
## See spec(...) for full column specifications.
#Read in Freedom Data
freedom <- read_excel("Freedom_In_The_50_States_2018.xlsx", sheet = "Overall")

#Combine
census_edu_2015 <- left_join(census_2015, edu_data_2015, by = c("CensusTract" = "Id2"))
full_2015 <- full_join(census_edu_2015, freedom %>% filter(Year == 2015), by = "State") 

#For the education data, we will need to divide the count by the total
full_2015_eduProp <- full_2015 %>% mutate(PropHighSchool = EstHighSchool/EstTotal,
                                          PropBachelors = EstBachelors/EstTotal,
                                          PropDoctorate = EstDoctorate/EstTotal)
full_2015_varOfInterest <- full_2015 %>%  select(IncomePerCap,Hispanic:Asian,Professional:Production, SelfEmployed, Unemployment, FiscalPolicy,EconomicFreedom, RegulatoryPolicy, PersonalFreedom, OverallFreedom)

“First Attempts” Geo Plots

full_2015_2 <- subset(full_2015, INTPTLONG > -150 & INTPTLONG < -50 & INTPTLAT > 24 & INTPTLAT < 50) 


loadPkg("ggplot2")
ggplot(full_2015_2, aes( x= INTPTLONG, y = INTPTLAT, color = Hispanic)) + 
  geom_point(size=0.005) + scale_color_gradient(low="orange", high="red") +
  labs (title="Scatter Plot: Spanish Ethnicity in 2015",
        x=" Latitude",
        y = "Longitude")

ggplot(full_2015_eduProp, aes( x= Professional, y = IncomePerCap, color = State == "New York", size = TotalPop)) +
  geom_point() +
  labs (title="Scatter Plot: IncomePerCap and Work type: Professional",
        x="Professional (% of pop in census tract)",
        y = "IncomePerCap")

# Population Histogram and QQ

loadPkg("ggplot2")

# NA's
dat0<-full_2015[!is.na(full_2015$IncomePerCap),]
sum(is.na(dat0$IncomePerCap))
## [1] 0
#removed outliers 
dat1 <- outlierKD2(dat0, IncomePerCap, TRUE)

## Outliers identified: 3589 
## Propotion (%) of outliers: 5.2 
## Mean of the outliers: 74140.26 
## Mean without removing outliers: 28491.23 
## Mean if we remove outliers: 26139.73 
## Outliers successfully removed
#Histogram of dependent variable
hist(dat1$TotalPop, 
        col = "#a8c2fb",
        main = "population Count",
        xlab = "US Dollars")

#ggplot
qqnorm(dat1$TotalPop, main="Q-Q plot of TotalPop")
qqline(dat1$TotalPop)

#Mean after removing NA's
format(mean(dat1$TotalPop, na.rm = TRUE))
## [1] "4369.254"
#Summary after removing NA's
summary(dat1$TotalPop)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       5    2929    4086    4369    5460   53812
#standard deviation
sd(dat1$TotalPop, na.rm = TRUE)
## [1] 2095.011

Income Histogram and QQ

#Histogram of dependent variable
hist(dat0$IncomePerCap, 
        col = "#a8c2fb",
        main = "Income per Capita",
        xlab = "US Dollars")

#ggplot
qqnorm(dat0$IncomePerCap, main="Q-Q plot of IncomePerCap")
qqline(dat0$IncomePerCap)

#Mean after removing NA's
format(mean(dat0$IncomePerCap, na.rm = TRUE))
## [1] "28491.23"
#Summary after removing NA's
summary(dat0$IncomePerCap)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     128   19123   25344   28491   33894  254204
#standard deviation
sd(dat0$IncomePerCap, na.rm = TRUE)
## [1] 15047.07
#Histogram of dependent variable
hist(dat1$IncomePerCap, 
        col = "#a8c2fb",
        main = "Income per Capita",
        xlab = "US Dollars")

#ggplot
qqnorm(dat1$IncomePerCap, main="Q-Q plot of IncomePerCap")
qqline(dat1$IncomePerCap)

#Mean after removing NA's
format(mean(dat1$IncomePerCap, na.rm = TRUE))
## [1] "26139.73"
#Summary after removing NA's
summary(dat1$IncomePerCap)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     128   18776   24730   26140   32247   56040    3589
#standard deviation
sd(dat1$IncomePerCap, na.rm = TRUE)
## [1] 10274.98

EDA By Group

loadPkg("dvmisc")
#Group by Economic Freedom
GroupEF <- quant_groups(dat1$EconomicFreedom, 4)
## Observations per group: 18417, 17756, 22780, 13244. 1064 missing.
str(GroupEF)
##  Factor w/ 4 levels "[-0.827,-0.256]",..: 3 3 3 3 3 3 3 3 3 3 ...
#Boxplot of IncomePerCap by Economic Freedom
plot(IncomePerCap ~ GroupEF, data=dat1, main="IncomePerCap and GroupEF", col=c("#ffd18b","#ffb76d","#ffa264", "#ff875f") )

summary(dat1$EconomicFreedom)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -0.8272 -0.2556  0.0152 -0.0866  0.1376  0.3550    1064
#Histogram of Economic Freedom
hist(dat1$EconomicFreedom, 
        col = "#ffd18b",
        main = "Economic Freedom")

#qqplot of Econoimc Freedom
qqnorm(dat1$EconomicFreedom, main="Q-Q plot of Economic Freedom")
qqline(dat1$EconomicFreedom)

#Group by PersonalFreedom
GroupPF <- quant_groups(dat1$PersonalFreedom, 4)
## Observations per group: 18061, 18704, 20888, 14544. 1064 missing.
str(GroupPF)
##  Factor w/ 4 levels "[-0.0444,0.0135]",..: 1 1 1 1 1 1 1 1 1 1 ...
#Boxplot of IncomePerCap by Personal Freedom
plot(IncomePerCap ~ GroupPF, data=dat1, main="IncomePerCap and GroupPF", col=c("#BFBFFF","#A3A3FF", "#7879FF", "#4949FF") )

summary(dat1$PersonalFreedom)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -0.0444  0.0135  0.0803  0.0641  0.1064  0.2450    1064
#Histogram of Personal Freedom
hist(dat1$PersonalFreedom, 
        col = "#A3A3FF",
        main = "Personal Freedom")

#qqplot of Personal Freedom
qqnorm(dat1$PersonalFreedom, main="Q-Q plot of Personal Freedom")
qqline(dat1$PersonalFreedom)

#Group by Regulatory Policy
GroupRP <- quant_groups(dat1$RegulatoryPolicy, 4)
## Observations per group: 19180, 17739, 17836, 17442. 1064 missing.
str(GroupRP)
##  Factor w/ 4 levels "[-0.457,-0.223]",..: 3 3 3 3 3 3 3 3 3 3 ...
#Boxplot of IncomePerCap by Regulatory Policy
plot(IncomePerCap ~ GroupRP, data=dat1, main="IncomePerCap and GroupRP", col=c("#FFDF01","#FED901", "#FFCF00", "#FEC300") )

summary(dat1$RegulatoryPolicy)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -0.4569 -0.2228 -0.0737 -0.1511 -0.0322  0.0715    1064
#Histogram of Regulatory Policy
hist(dat1$RegulatoryPolicy, 
        col = "#FFDF01",
        main = "Regulatory Policy ")

#qqplot of Regulatory Policy
qqnorm(dat1$RegulatoryPolicy, main="Q-Q plot of Regulatory Policy")
qqline(dat1$RegulatoryPolicy)

#Group by Fiscal Policy
GroupFP <- quant_groups(dat1$FiscalPolicy, 4)
## Observations per group: 18164, 19958, 16580, 17495. 1064 missing.
str(GroupFP)
##  Factor w/ 4 levels "[-0.37,-0.0602]",..: 3 3 3 3 3 3 3 3 3 3 ...
#Boxplot of IncomePerCap by Fiscal Policy
plot(IncomePerCap ~ GroupFP, data=dat1, main="IncomePerCap and GroupFP", col=c("#7EFFD4","#70EBBA", "#64D8A7", "#58BD95") )

summary(dat1$FiscalPolicy)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -0.3702 -0.0602  0.0634  0.0646  0.1767  0.4024    1064
#Histogram of Fiscal Policy
hist(dat1$FiscalPolicy, 
        col = "#7EFFD4",
        main = "Fiscal Policy")

#qqplot of Fiscal Policy
qqnorm(dat1$FiscalPolicy, main="Q-Q plot of Fiscal Policy")
qqline(dat1$FiscalPolicy)

#Group by Overall Freedom
GroupOF <- quant_groups(dat1$OverallFreedom, 4)
## Observations per group: 18440, 17906, 18395, 17456. 1064 missing.
str(GroupOF)
##  Factor w/ 4 levels "[-0.814,-0.102]",..: 2 2 2 2 2 2 2 2 2 2 ...
#Boxplot of IncomePerCap by Overall Freedom
plot(IncomePerCap ~ GroupOF, data=dat1, main="IncomePerCap and GroupOF", col=c("#C0C6CB","#98A5C0", "#7688BB", "#536CB5") )

summary(dat1$OverallFreedom)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
## -0.8136 -0.1021  0.0652 -0.0224  0.1637  0.4614    1064
#Histogram of Overall Freedom
hist(dat1$OverallFreedom, 
        col = "#A3A3FF",
        main = "Overall Freedom")

#qqplot of Overall Freedom
qqnorm(dat1$OverallFreedom, main="Q-Q plot of Overall Freedom")
qqline(dat1$OverallFreedom)

#Group by Unemployment
GroupUnemployment <- quant_groups(dat1$Unemployment, 4)
## Observations per group: 18765, 18330, 18095, 17970. 101 missing.
str(GroupUnemployment)
##  Factor w/ 4 levels "[0,5.1]","(5.1,7.7]",..: 2 4 2 3 1 3 3 3 3 2 ...
#Boxplot of IncomePerCap by Unemployment
plot(IncomePerCap ~ GroupUnemployment, data=dat1, main="IncomePerCap and GroupUnemployment", col=c("#FFF57B","#FFE469", "#FECC51", "#FCB033") )

summary(dat1$Unemployment)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   5.100   7.700   9.028  11.400 100.000     101
#Histogram of Unemployment
hist(dat1$Unemployment, 
        col = "#A3A3FF",
       )

#qqplot of Unemployment Unemployment
qqnorm(dat1$Unemployment, main="Q-Q plot of Unemployment")
qqline(dat1$Unemployment)

#Group by Professional
GroupProfessional <- quant_groups(dat1$Professional, 4)
## Observations per group: 18379, 18323, 18173, 18281. 105 missing.
str(GroupProfessional)
##  Factor w/ 4 levels "[0,24.1]","(24.1,32.6]",..: 3 1 2 2 4 2 1 3 2 2 ...
#Boxplot of IncomePerCap by Professional
plot(IncomePerCap ~ GroupProfessional, data=dat1, main="IncomePerCap and GroupProfessional", col=c("#F3D8F2","#E6B2E4", "#D88DD5", "#CA68C7") )

summary(dat1$Professional)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0    24.1    32.6    34.8    43.8   100.0     105
#Histogram of Professional
hist(dat1$Professional, 
        col = "#D88DD5",
        main = "Professional")

#qqplot of Professional
qqnorm(dat1$Professional, main="Q-Q plot of Professional")
qqline(dat1$Professional)

#Group by Office
GroupOffice <- quant_groups(dat1$Office, 4)
## Observations per group: 18303, 18828, 17766, 18259. 105 missing.
str(GroupOffice)
##  Factor w/ 4 levels "[0,20.1]","(20.1,23.8]",..: 2 2 2 3 1 4 3 4 3 1 ...
#Boxplot of IncomePerCap by Office
plot(IncomePerCap ~ GroupOffice, data=dat1, main="IncomePerCap and GroupOffice", col=c("#D1FFD5","#B4FFB2", "#98FF98", "#79F58A") )

summary(dat1$Office)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   20.10   23.80   23.95   27.50  100.00     105
#Histogram of Office
hist(dat1$Office, 
        col = "#FECC51",
        main = "Office")

#qqplot of Office
qqnorm(dat1$Office, main="Q-Q plot of Office")
qqline(dat1$Office)

#Group by Service
GroupService <- quant_groups(dat1$Service, 4)
## Observations per group: 18672, 18068, 18275, 18141. 105 missing.
str(GroupService)
##  Factor w/ 4 levels "[0,13.5]","(13.5,17.9]",..: 2 4 4 3 2 2 4 1 2 2 ...
#Boxplot of IncomePerCap by Service
plot(IncomePerCap ~ GroupService, data=dat1, main="IncomePerCap and GroupService", col=c("#E0BCBF","#D8ABB1", "#CF989F", "#C0838C") )

summary(dat1$Service)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0    13.5    17.9    19.1    23.6   100.0     105
#Histogram of Service
hist(dat1$Service, 
        col = "#79F58A",
        main = "Service")

#qqplot of Service
qqnorm(dat1$Service, main="Q-Q plot of Service")
qqline(dat1$Service)

#Group by Construction
GroupConstruction <- quant_groups(dat1$Construction, 4)
## Observations per group: 18588, 18209, 18113, 18246. 105 missing.
str(GroupConstruction)
##  Factor w/ 4 levels "[0,5]","(5,8.4]",..: 3 3 3 3 1 2 3 2 2 3 ...
#Boxplot of IncomePerCap by Construction
plot(IncomePerCap ~ GroupConstruction, data=dat1, main="IncomePerCap and GroupConstruction", col=c("#A9AB98","#949180", "#7D7968", "#5E594F") )

summary(dat1$Construction)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   5.000   8.400   9.295  12.500 100.000     105
#Histogram of Construciton
hist(dat1$Construction, 
        col = "#C0838C",
        main = "Construction")

#qqplot of Construction
qqnorm(dat1$Construction, main="Q-Q plot of Construction")
qqline(dat1$Construction)

#Group by Production
GroupProduction <- quant_groups(dat1$Production, 4)
## Observations per group: 18566, 18226, 18085, 18279. 105 missing.
str(GroupProduction)
##  Factor w/ 4 levels "[0,7.1]","(7.1,11.8]",..: 3 4 3 3 3 3 3 1 3 4 ...
#Boxplot of IncomePerCap by Production
plot(IncomePerCap ~ GroupProduction, data=dat1, main="IncomePerCap and Production", col=c("#F3D8F2","#E6B2E4", "#D88DD5", "#CA68C7") )

summary(dat1$Production)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00    7.10   11.80   12.86   17.40  100.00     105
#Histogram of Producrtion
hist(dat1$Production, 
        col = "#E6B2E4",
        main = "Production")

#qqplot of Production
qqnorm(dat1$Production, main="Q-Q plot of Production")
qqline(dat1$Production)

#Group by Self-Employed
GroupProduction <- quant_groups(dat1$SelfEmployed, 4)
## Observations per group: 19155, 17839, 18195, 17967. 105 missing.
str(GroupProduction)
##  Factor w/ 4 levels "[0,3.6]","(3.6,5.5]",..: 2 3 4 1 2 3 1 3 2 3 ...
#Boxplot of IncomePerCap by Selfemployed
plot(IncomePerCap ~ GroupProduction, data=dat1, main="IncomePerCap and Selfemployed", col=c("#F3D8F2","#E6B2E4", "#D88DD5", "#CA68C7") )

summary(dat1$SelfEmployed)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   3.600   5.500   6.227   8.100 100.000     105
#Histogram of Self Employed
hist(dat1$SelfEmployed, 
        col = "#E6B2E4",
        main = "Production")

#qqplot of Self employed
qqnorm(dat1$SelfEmployed, main="Q-Q plot of Selfemployed")
qqline(dat1$SelfEmployed)

#Group by Black
GroupBlack <- quant_groups(dat1$Black, 4)
## Observations per group: 18430, 18252, 18303, 18276. 0 missing.
str(GroupBlack)
##  Factor w/ 4 levels "[0,0.7]","(0.7,3.7]",..: 3 4 4 2 4 3 4 3 3 3 ...
#Boxplot of IncomePerCap by Black Population
plot(IncomePerCap ~ GroupBlack, data=dat1, main="IncomePerCap and GroupBlack", col=c("#d0dfff","#a8c2fb", "#86abf9", "#6893ee") )

summary(dat1$Black)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.70    3.70   13.27   14.40  100.00
#Histogram of Black
hist(dat1$Black, 
        col = "#d0dfff",
        main = "Black")

#qqplot of Black
qqnorm(dat1$Black, main="Q-Q plot of Black")
qqline(dat1$Black)

#Group by Hispanic
GroupHispanic <- quant_groups(dat1$Hispanic, 4)
## Observations per group: 18472, 18246, 18233, 18310. 0 missing.
str(GroupHispanic)
##  Factor w/ 4 levels "[0,2.4]","(2.4,7]",..: 1 1 1 3 1 3 2 1 1 1 ...
#Boxplot of IncomePerCap by Hispanic Population
plot(IncomePerCap ~ GroupHispanic, data=dat1, main="IncomePerCap and GroupHispanic", col=c("#F6BDC0","#F1959B", "#F07470", "#EA4C46") )

summary(dat1$Hispanic)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    2.40    7.00   16.86   20.40  100.00
#Histogram of Hispanic
hist(dat1$Hispanic, 
        col = "#EA4C46",
        main = "Hispanic")

#qqplot of Hispanic
qqnorm(dat1$Hispanic, main="Q-Q plot of Hspanic")
qqline(dat1$Hispanic)

#Group by Asians
GroupAsian <- quant_groups(dat1$Asian, 3)
## Observations per group: 26078, 23036, 24147. 0 missing.
str(GroupAsian)
##  Factor w/ 3 levels "[0,0.5]","(0.5,3.2]",..: 2 2 2 1 3 1 1 1 1 1 ...
#Boxplot of IncomePerCap by Asian Population
plot(IncomePerCap ~ GroupAsian, data=dat1, main="IncomePerCap and GroupAsian", col=c("#FFE6C8","#FFCCBE", "#EFABA0", "#D6806F") )

summary(dat1$Asian)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.20    1.40    4.59    4.80   91.30
#Histogram of Asian
hist(dat1$Asian, 
        col = "#FFE6C8",
        main = "Asian")

#qqplot of Asian
qqnorm(dat1$Asian, main="Q-Q plot of Asian")
qqline(dat1$Asian)

#Group by White
GroupWhite <- quant_groups(dat1$White, 4)
## Observations per group: 18351, 18331, 18285, 18294. 0 missing.
str(GroupWhite)
##  Factor w/ 4 levels "[0,39.4]","(39.4,71.4]",..: 3 2 3 3 2 3 3 3 4 3 ...
#Boxplot of IncomePerCap by White Population
plot(IncomePerCap ~ GroupWhite, data=dat1, main="IncomePerCap and GroupWhite", col=c("#F5F8FA", "#EFF2F4", "#EAEDEF", "#E0E3E5") )

summary(dat1$White)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00   39.40   71.40   62.03   88.30  100.00
#Histogram of White
hist(dat1$White, 
        col = "#F5F8FA",
        main = "White")

#qqplot of White
qqnorm(dat1$White, main="Q-Q plot of White")
qqline(dat1$White)

Correlations

full_2015_cor <- cor(full_2015_varOfInterest, use = "complete.obs")

loadPkg("corrplot")
## corrplot 0.84 loaded
corrplot(full_2015_cor)

Freedom_2015 <- subset(full_2015_varOfInterest, select = c(IncomePerCap, OverallFreedom, EconomicFreedom, PersonalFreedom, RegulatoryPolicy, FiscalPolicy))

Work_2015 <- subset(full_2015_varOfInterest, select = c(IncomePerCap, Professional, Service, Office, Construction, Production, SelfEmployed, Unemployment))

Ethnic_2015 <- subset(full_2015_varOfInterest, select = c(IncomePerCap, Hispanic, White, Black, Native, Asian))
Freedom_2015_cor <- cor(Freedom_2015, use = "complete.obs")
#corrplot(Freedom_2015_cor)
corrplot.mixed(Freedom_2015_cor, tl.pos = "lt")

Work_2015_cor <- cor(Work_2015, use = "complete.obs")
#corrplot(Work_2015_cor)
corrplot.mixed(Work_2015_cor, tl.pos = "lt")

Ethnic_2015_cor <- cor(Ethnic_2015, use = "complete.obs")
#corrplot(Ethnic_2015_cor)
corrplot.mixed(Ethnic_2015_cor, tl.pos = "lt")

Scatter PLots

ggplot(dat1, aes( x= Professional, y = IncomePerCap, color = White)) + 
  geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~Professional) +
  labs (title="Scatter Plot of IncomePerCap and Professional",
        x="Professional",
        y = "IncomePerCap")

ggplot(dat1, aes( x= Unemployment, y = IncomePerCap, color = White)) + 
  geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~Unemployment) +
  labs (title="Scatter Plot of IncomePerCap and Unemployment",
        x="Unemployment",
        y = "IncomePerCap")

ggplot(dat1, aes( x= Service, y = IncomePerCap, color = White)) + 
  geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~Service) +
  labs (title="Scatter Plot of IncomePerCap and Service",
        x="Service",
        y = "IncomePerCap")

ggplot(dat1, aes( x= EconomicFreedom, y = IncomePerCap, color = White)) + 
  geom_point(size=0.1) + geom_smooth(method='lm', formula= IncomerPerCap~EconomicFreedom) +
  labs (title="Scatter Plot of IncomePerCap and EconomicFreedom",
        x="EconomicFreedom",
        y = "IncomePerCap")

ANOVA

#' @title Returns the max value of each row of a data.frame or matrix
#'
#' @description
#' Returns maximum value of each row of a data.frame or matrix.
#' @param df Data.frame or matrix, required.
#' @param na.rm Logical value, optional, TRUE by default. Defines whether NA values should be removed first. Otherwise result will be NA when any NA is in the given vector.
#' @return Returns a vector of numbers of length equal to number of rows in df.
#' @template maxmin
#' @export
rowMaxs <- function(df, na.rm=TRUE) {

  if (is.matrix(df)) {df <- data.frame(df, stringsAsFactors=FALSE, drop = FALSE)}

  valid.cols <- sapply(df, function(x) { is.numeric(x) || is.logical(x) || is.character(x)})
  stopifnot(any(valid.cols))
  # or could just return NA?:
  # if (!any(valid.cols)) {return(NA)}
  if (any(!valid.cols) ) {warning('using only numeric (double or integer) or logical or character columns -- ignoring other columns ')}

  result <- do.call(pmax, c(df[ , valid.cols, drop = FALSE], na.rm=na.rm))

  result[nononmissing <- rowSums(!is.na(df[ , valid.cols, drop = FALSE]))==0] <- -Inf
  if (any(nononmissing)) {warning('where no non-missing arguments, returning -Inf')}
  return(result)

  # df = data.frame of numeric values, i.e. a list of vectors passed to pmax
  # Value returned is vector, each element is max of a row of df
}

full_2015_varOfInterest1 <- full_2015_varOfInterest %>% mutate(EthnicMajority = case_when(
  White >= 50 ~ "White",
  Hispanic >= 50 ~ "Hispanic",
  Black >= 50 ~ "Black",
  Native >= 50 ~ "Native",
  Asian >= 50 ~ "Asian",
  TRUE ~ "Diverse"))
EthnicityMaxes <- rowMaxs(full_2015_varOfInterest1 %>%  select(White, Hispanic, Black, Native, Asian))
WorkMaxes <- rowMaxs(full_2015_varOfInterest1 %>%  select(Professional:Unemployment))
full_2015_varOfInterest2 <- cbind(full_2015_varOfInterest1, EthnicityMaxes, WorkMaxes)

full_2015_varOfInterest_withMajorityEthnicity <- full_2015_varOfInterest2 %>% mutate(EthnicPlurality = case_when(
  White == EthnicityMaxes ~ "White",
  Hispanic == EthnicityMaxes ~ "Hispanic",
  Black == EthnicityMaxes ~ "Black",
  Native == EthnicityMaxes ~ "Native",
  Asian == EthnicityMaxes ~ "Asian",
  TRUE ~ "Error"),
  WorkPlurality = case_when(
  Professional == WorkMaxes ~ "Professional",
  Service == WorkMaxes ~ "Service",
  Office == WorkMaxes ~ "Office",
  Construction == WorkMaxes ~ "Construction",
  Production == WorkMaxes ~ "Production",
  SelfEmployed == WorkMaxes ~ "SelfEmployed",
  Unemployment == WorkMaxes ~ "Unemployment",
  TRUE ~ "Error"))
full_2015_varOfInterest_withMajorityEthnicity$EthnicMajority <-
  as.factor(full_2015_varOfInterest_withMajorityEthnicity$EthnicMajority)
full_2015_varOfInterest_withMajorityEthnicity$EthnicPlurality <-
  as.factor(full_2015_varOfInterest_withMajorityEthnicity$EthnicPlurality)
full_2015_varOfInterest_withMajorityEthnicity$WorkPlurality <-
  as.factor(full_2015_varOfInterest_withMajorityEthnicity$WorkPlurality)
anova_dat <- full_2015_varOfInterest_withMajorityEthnicity %>% filter(EthnicPlurality!= "Error",
                                                                 WorkPlurality!= "Error")
Ethnic_2015_anova = aov(IncomePerCap ~ EthnicPlurality, data = anova_dat)
Ethnic_2015_anova
## Call:
##    aov(formula = IncomePerCap ~ EthnicPlurality, data = anova_dat)
## 
## Terms:
##                 EthnicPlurality    Residuals
## Sum of Squares     2.590039e+12 1.394545e+13
## Deg. of Freedom               4        73155
## 
## Residual standard error: 13806.84
## Estimated effects may be unbalanced
## 39 observations deleted due to missingness
names(Ethnic_2015_anova)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "na.action"     "contrasts"     "xlevels"       "call"         
## [13] "terms"         "model"
summary(Ethnic_2015_anova)
##                    Df    Sum Sq   Mean Sq F value Pr(>F)    
## EthnicPlurality     4 2.590e+12 6.475e+11    3397 <2e-16 ***
## Residuals       73155 1.395e+13 1.906e+08                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 39 observations deleted due to missingness
ggplot(anova_dat, aes(x=EthnicPlurality, y=IncomePerCap)) + 
  geom_boxplot(outlier.shape=8, outlier.size=4) +
  labs(title="Income/Capita with Different Ethnic Majority", x="Ethnic Majority", y = "Income per Cap")

TukeyHSD(Ethnic_2015_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = IncomePerCap ~ EthnicPlurality, data = anova_dat)
## 
## $EthnicPlurality
##                        diff         lwr          upr     p adj
## Black-Asian     -12780.4809 -13907.7147 -11653.24708 0.0000000
## Hispanic-Asian  -13705.6676 -14812.6303 -12598.70484 0.0000000
## Native-Asian    -16464.3603 -19310.6485 -13618.07205 0.0000000
## White-Asian        647.5161   -403.9082   1698.94040 0.4464321
## Hispanic-Black    -925.1867  -1505.7776   -344.59577 0.0001346
## Native-Black     -3683.8794  -6369.5964   -998.16235 0.0017134
## White-Black      13427.9970  12961.9366  13894.05745 0.0000000
## Native-Hispanic  -2758.6927  -5435.9649    -81.42052 0.0396499
## White-Hispanic   14353.1837  13938.5480  14767.81938 0.0000000
## White-Native     17111.8764  14457.0858  19766.66697 0.0000000
Work_2015_anova = aov(IncomePerCap ~ WorkPlurality, data = anova_dat)
Work_2015_anova
## Call:
##    aov(formula = IncomePerCap ~ WorkPlurality, data = anova_dat)
## 
## Terms:
##                 WorkPlurality    Residuals
## Sum of Squares   4.242548e+12 1.229294e+13
## Deg. of Freedom             6        73153
## 
## Residual standard error: 12963.19
## Estimated effects may be unbalanced
## 39 observations deleted due to missingness
names(Work_2015_anova)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "na.action"     "contrasts"     "xlevels"       "call"         
## [13] "terms"         "model"
summary(Work_2015_anova)
##                  Df    Sum Sq   Mean Sq F value Pr(>F)    
## WorkPlurality     6 4.243e+12 7.071e+11    4208 <2e-16 ***
## Residuals     73153 1.229e+13 1.680e+08                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 39 observations deleted due to missingness
ggplot(anova_dat, aes(x=WorkPlurality, y=IncomePerCap)) + 
  geom_boxplot(outlier.shape=8, outlier.size=4) +
  labs(title="Income/Capita with Different Work Majority", x="Work Majority", y = "Income per Cap")

TukeyHSD(Work_2015_anova)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = IncomePerCap ~ WorkPlurality, data = anova_dat)
## 
## $WorkPlurality
##                                   diff           lwr         upr     p adj
## Office-Construction         4152.88357   2905.526990   5400.2402 0.0000000
## Production-Construction     1332.84270     -5.891578   2671.5770 0.0518977
## Professional-Construction  17590.26271  16385.932493  18794.5929 0.0000000
## SelfEmployed-Construction  11446.99976    349.799288  22544.2002 0.0380499
## Service-Construction          73.21197  -1179.886893   1326.3108 0.9999978
## Unemployment-Construction  -5476.10244  -7681.711665  -3270.4932 0.0000000
## Production-Office          -2820.04087  -3533.458145  -2106.6236 0.0000000
## Professional-Office        13437.37913  13028.519743  13846.2385 0.0000000
## SelfEmployed-Office         7294.11618  -3745.114449  18333.3468 0.4482729
## Service-Office             -4079.67161  -4615.406142  -3543.9371 0.0000000
## Unemployment-Office        -9628.98602 -11521.462377  -7736.5097 0.0000000
## Professional-Production    16257.42001  15622.221597  16892.6184 0.0000000
## SelfEmployed-Production    10114.15706   -935.771630  21164.0857 0.0985074
## Service-Production         -1259.63073  -1983.041067   -536.2204 0.0000059
## Unemployment-Production    -6808.94514  -8762.858599  -4855.0317 0.0000000
## SelfEmployed-Professional  -6143.26295 -17177.714716   4891.1888 0.6552473
## Service-Professional      -17517.05074 -17943.107437 -17090.9940 0.0000000
## Unemployment-Professional -23066.36515 -24930.763066 -21201.9672 0.0000000
## Service-SelfEmployed      -11373.78779 -22413.668735   -333.9068 0.0384811
## Unemployment-SelfEmployed -16923.10220 -28111.195272  -5735.0091 0.0001665
## Unemployment-Service       -5549.31441  -7445.580501  -3653.0483 0.0000000

Chi Squared Tests

loadPkg('sjPlot')
## Learn more about sjPlot with 'browseVignettes("sjPlot")'.
loadPkg('visualize')

chi_table_dat <- full_2015_varOfInterest_withMajorityEthnicity %>%
  select(c(WorkPlurality, EthnicPlurality)) %>% filter(WorkPlurality != "Error")

chi_table_dat %>%
  sjtab(fun = "xtab", var.labels=c("Work", "Ethnicity"),
       show.summary=T, show.exp=T, show.legend=T)
Work Ethnicity Total
Asian Black Error Hispanic Native White
Construction 0
19
31
104
0
0
639
137
1
3
359
767
1030
1030
Error 0
0
0
0
0
0
0
0
0
0
0
0
0
0
Office 128
193
1564
1087
0
0
2199
1426
17
30
6816
7989
10724
10724
Production 19
70
480
398
0
0
993
522
6
11
2426
2923
3924
3924
Professional 927
851
2346
4804
0
0
2510
6299
114
131
41484
35296
47381
47381
SelfEmployed 0
0
0
1
0
0
2
2
0
0
11
10
13
13
Service 241
174
2744
984
0
0
3274
1290
49
27
3394
7227
9702
9702
Unemployment 0
8
257
43
0
0
114
56
15
1
39
317
425
425
Total 1315
1315
7422
7422
0
0
9731
9731
202
202
54529
54529
73199
73199
χ2=NaN · df=35 · Cramer’s V=NaN · Fisher’s p=0.000

observed values
expected values

chi.test <- chisq.test(full_2015_varOfInterest_withMajorityEthnicity$EthnicPlurality, full_2015_varOfInterest_withMajorityEthnicity$WorkPlurality)

#Exhaustive Search

loadPkg('leaps')
#This is essentially best fit 
reg.best10 <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17)  
## Reordering variables and trying again:
plot(reg.best10, scale = "adjr2", main = "Adjusted R^2")

plot(reg.best10, scale = "r2", main = "R^2")

# In the "leaps" package, we can use scale=c("bic","Cp","adjr2","r2")
plot(reg.best10, scale = "bic", main = "BIC")

plot(reg.best10, scale = "Cp", main = "Cp")

reg.summary<-summary(reg.best10)
names(reg.best10)
##  [1] "np"        "nrbar"     "d"         "rbar"      "thetab"    "first"    
##  [7] "last"      "vorder"    "tol"       "rss"       "bound"     "nvmax"    
## [13] "ress"      "ir"        "nbest"     "lopt"      "il"        "ier"      
## [19] "xnames"    "method"    "force.in"  "force.out" "sserr"     "intercept"
## [25] "lindep"    "reorder"   "nullrss"   "nn"        "call"
plot(reg.summary$rsq,xlab = "# of Varibales", ylab = "Rsquare", type = "l")

plot(reg.summary$rss,xlab = "# of Varibales", ylab = "RSS", type = "l")

#Adjusted R2
which.max(reg.summary$adjr2)
## [1] 14
plot(reg.summary$adjr2,xlab = '# of Varibales', ylab = 'Adjusted Rsq', type = "l")
points(12,reg.summary$adjr2[12],col = "red",cex = 2, pch =20)

#CP
which.min(reg.summary$cp)
## [1] 13
plot(reg.summary$cp,xlab = '# of Varibales', ylab = 'CP', type = "l")
points(12,reg.summary$cp[12],col = "red",cex = 2, pch =20)

#BIC
which.min(reg.summary$bic)
## [1] 12
plot(reg.summary$bic,xlab = '# of Varibales', ylab = 'BIC', type = "l")
points(10,reg.summary$bic[12],col = "red",cex = 2, pch =20)

Forward Selection

reg.forward10 <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17, nbest = 2, method = "forward")
## Reordering variables and trying again:
plot(reg.forward10, scale = "adjr2", main = "Adjusted R^2")

plot(reg.forward10, scale = "bic", main = "BIC")

plot(reg.forward10, scale = "Cp", main = "Cp")

#summary(reg.forward10)

Backward Selection

Now backwards (nvmax=10 and nbest=2)

reg.back10 <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17, nbest = 2, method = "backward")
## Reordering variables and trying again:
plot(reg.back10, scale = "adjr2", main = "Adjusted R^2")

plot(reg.back10, scale = "bic", main = "BIC")

plot(reg.back10, scale = "Cp", main = "Cp")

summary(reg.back10)
## Subset selection object
## Call: regsubsets.formula(IncomePerCap ~ ., data = full_2015_varOfInterest, 
##     nvmax = 17, nbest = 2, method = "backward")
## 17 Variables  (and intercept)
##                  Forced in Forced out
## Hispanic             FALSE      FALSE
## White                FALSE      FALSE
## Black                FALSE      FALSE
## Native               FALSE      FALSE
## Asian                FALSE      FALSE
## Professional         FALSE      FALSE
## Service              FALSE      FALSE
## Office               FALSE      FALSE
## Construction         FALSE      FALSE
## Production           FALSE      FALSE
## SelfEmployed         FALSE      FALSE
## Unemployment         FALSE      FALSE
## FiscalPolicy         FALSE      FALSE
## EconomicFreedom      FALSE      FALSE
## PersonalFreedom      FALSE      FALSE
## RegulatoryPolicy     FALSE      FALSE
## OverallFreedom       FALSE      FALSE
## 2 subsets of each size up to 15
## Selection Algorithm: backward
##           Hispanic White Black Native Asian Professional Service Office
## 1  ( 1 )  " "      " "   " "   " "    " "   "*"          " "     " "   
## 1  ( 2 )  " "      "*"   " "   " "    " "   " "          " "     " "   
## 2  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 2  ( 2 )  " "      "*"   "*"   " "    " "   " "          " "     " "   
## 3  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 3  ( 2 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 4  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 4  ( 2 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 5  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 5  ( 2 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 6  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 6  ( 2 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 7  ( 1 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 7  ( 2 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 8  ( 1 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 8  ( 2 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 9  ( 1 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 9  ( 2 )  " "      "*"   "*"   " "    " "   "*"          "*"     " "   
## 10  ( 1 ) " "      "*"   "*"   " "    " "   "*"          "*"     " "   
## 10  ( 2 ) "*"      "*"   "*"   " "    " "   "*"          "*"     " "   
## 11  ( 1 ) "*"      "*"   "*"   " "    " "   "*"          "*"     " "   
## 11  ( 2 ) "*"      "*"   "*"   "*"    " "   "*"          "*"     " "   
## 12  ( 1 ) "*"      "*"   "*"   "*"    " "   "*"          "*"     " "   
## 12  ( 2 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     " "   
## 13  ( 1 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     " "   
## 13  ( 2 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     "*"   
## 14  ( 1 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     " "   
## 14  ( 2 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     "*"   
## 15  ( 1 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     "*"   
##           Construction Production SelfEmployed Unemployment FiscalPolicy
## 1  ( 1 )  " "          " "        " "          " "          " "         
## 1  ( 2 )  " "          " "        " "          " "          " "         
## 2  ( 1 )  " "          " "        " "          " "          " "         
## 2  ( 2 )  " "          " "        " "          " "          " "         
## 3  ( 1 )  " "          " "        "*"          " "          " "         
## 3  ( 2 )  " "          " "        " "          " "          " "         
## 4  ( 1 )  " "          " "        " "          " "          "*"         
## 4  ( 2 )  " "          " "        "*"          "*"          " "         
## 5  ( 1 )  " "          " "        " "          "*"          "*"         
## 5  ( 2 )  " "          " "        "*"          "*"          " "         
## 6  ( 1 )  " "          " "        "*"          "*"          "*"         
## 6  ( 2 )  " "          " "        "*"          "*"          "*"         
## 7  ( 1 )  " "          " "        "*"          "*"          "*"         
## 7  ( 2 )  "*"          "*"        "*"          "*"          " "         
## 8  ( 1 )  "*"          " "        "*"          "*"          "*"         
## 8  ( 2 )  "*"          "*"        "*"          "*"          "*"         
## 9  ( 1 )  "*"          "*"        "*"          "*"          "*"         
## 9  ( 2 )  "*"          "*"        "*"          "*"          "*"         
## 10  ( 1 ) "*"          "*"        "*"          "*"          "*"         
## 10  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 11  ( 1 ) "*"          "*"        "*"          "*"          "*"         
## 11  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 12  ( 1 ) "*"          "*"        "*"          "*"          "*"         
## 12  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 13  ( 1 ) "*"          "*"        "*"          "*"          "*"         
## 13  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 14  ( 1 ) "*"          "*"        "*"          "*"          "*"         
## 14  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 15  ( 1 ) "*"          "*"        "*"          "*"          "*"         
##           EconomicFreedom RegulatoryPolicy PersonalFreedom OverallFreedom
## 1  ( 1 )  " "             " "              " "             " "           
## 1  ( 2 )  " "             " "              " "             " "           
## 2  ( 1 )  " "             " "              " "             " "           
## 2  ( 2 )  " "             " "              " "             " "           
## 3  ( 1 )  " "             " "              " "             " "           
## 3  ( 2 )  "*"             " "              " "             " "           
## 4  ( 1 )  "*"             " "              " "             " "           
## 4  ( 2 )  " "             " "              " "             " "           
## 5  ( 1 )  "*"             " "              " "             " "           
## 5  ( 2 )  " "             " "              " "             " "           
## 6  ( 1 )  "*"             " "              " "             " "           
## 6  ( 2 )  " "             " "              " "             " "           
## 7  ( 1 )  "*"             " "              " "             " "           
## 7  ( 2 )  " "             " "              " "             " "           
## 8  ( 1 )  "*"             " "              " "             " "           
## 8  ( 2 )  " "             " "              " "             " "           
## 9  ( 1 )  "*"             " "              " "             " "           
## 9  ( 2 )  " "             " "              " "             " "           
## 10  ( 1 ) "*"             " "              " "             " "           
## 10  ( 2 ) " "             " "              " "             " "           
## 11  ( 1 ) "*"             " "              " "             " "           
## 11  ( 2 ) " "             " "              " "             " "           
## 12  ( 1 ) "*"             " "              " "             " "           
## 12  ( 2 ) " "             " "              " "             " "           
## 13  ( 1 ) "*"             " "              " "             " "           
## 13  ( 2 ) " "             " "              " "             " "           
## 14  ( 1 ) "*"             " "              "*"             " "           
## 14  ( 2 ) "*"             " "              " "             " "           
## 15  ( 1 ) "*"             " "              "*"             " "

Sequential Replacement seqrep

reg.seqrep <- regsubsets(IncomePerCap~. , data = full_2015_varOfInterest, nvmax = 17, nbest = 2 , method = "seqrep")
## Reordering variables and trying again:
plot(reg.seqrep, scale = "adjr2", main = "Adjusted R^2")

plot(reg.seqrep, scale = "bic", main = "BIC")

plot(reg.seqrep, scale = "Cp", main = "Cp")

summary(reg.seqrep)
## Subset selection object
## Call: regsubsets.formula(IncomePerCap ~ ., data = full_2015_varOfInterest, 
##     nvmax = 17, nbest = 2, method = "seqrep")
## 17 Variables  (and intercept)
##                  Forced in Forced out
## Hispanic             FALSE      FALSE
## White                FALSE      FALSE
## Black                FALSE      FALSE
## Native               FALSE      FALSE
## Asian                FALSE      FALSE
## Professional         FALSE      FALSE
## Service              FALSE      FALSE
## Office               FALSE      FALSE
## Construction         FALSE      FALSE
## Production           FALSE      FALSE
## SelfEmployed         FALSE      FALSE
## Unemployment         FALSE      FALSE
## FiscalPolicy         FALSE      FALSE
## EconomicFreedom      FALSE      FALSE
## PersonalFreedom      FALSE      FALSE
## RegulatoryPolicy     FALSE      FALSE
## OverallFreedom       FALSE      FALSE
## 2 subsets of each size up to 15
## Selection Algorithm: 'sequential replacement'
##           Hispanic White Black Native Asian Professional Service Office
## 1  ( 1 )  " "      " "   " "   " "    " "   "*"          " "     " "   
## 1  ( 2 )  " "      " "   " "   " "    " "   " "          "*"     " "   
## 2  ( 1 )  " "      " "   " "   " "    " "   "*"          " "     " "   
## 2  ( 2 )  " "      " "   "*"   " "    " "   "*"          " "     " "   
## 3  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 3  ( 2 )  " "      " "   " "   " "    " "   "*"          " "     " "   
## 4  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 4  ( 2 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 5  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     " "   
## 5  ( 2 )  " "      " "   " "   " "    " "   "*"          "*"     " "   
## 6  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     "*"   
## 6  ( 2 )  " "      "*"   " "   " "    " "   "*"          "*"     " "   
## 7  ( 1 )  " "      "*"   " "   " "    " "   "*"          " "     "*"   
## 7  ( 2 )  " "      "*"   " "   " "    " "   "*"          " "     "*"   
## 8  ( 1 )  " "      "*"   " "   " "    " "   "*"          "*"     "*"   
## 8  ( 2 )  " "      "*"   " "   " "    " "   "*"          "*"     "*"   
## 9  ( 1 )  " "      "*"   " "   "*"    " "   "*"          "*"     "*"   
## 9  ( 2 )  " "      "*"   " "   "*"    " "   "*"          "*"     "*"   
## 10  ( 1 ) " "      "*"   " "   "*"    "*"   "*"          "*"     "*"   
## 10  ( 2 ) " "      "*"   " "   "*"    "*"   "*"          "*"     "*"   
## 11  ( 1 ) " "      "*"   " "   "*"    "*"   "*"          "*"     " "   
## 11  ( 2 ) " "      "*"   " "   "*"    "*"   "*"          "*"     " "   
## 12  ( 1 ) "*"      "*"   "*"   "*"    " "   "*"          "*"     " "   
## 12  ( 2 ) "*"      "*"   "*"   "*"    " "   "*"          "*"     " "   
## 13  ( 1 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     " "   
## 13  ( 2 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     " "   
## 14  ( 1 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     " "   
## 14  ( 2 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     " "   
## 15  ( 1 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     "*"   
## 15  ( 2 ) "*"      "*"   "*"   "*"    "*"   "*"          "*"     "*"   
##           Construction Production SelfEmployed Unemployment FiscalPolicy
## 1  ( 1 )  " "          " "        " "          " "          " "         
## 1  ( 2 )  " "          " "        " "          " "          " "         
## 2  ( 1 )  " "          " "        " "          "*"          " "         
## 2  ( 2 )  " "          " "        " "          " "          " "         
## 3  ( 1 )  " "          " "        " "          " "          " "         
## 3  ( 2 )  " "          " "        " "          "*"          " "         
## 4  ( 1 )  " "          " "        " "          "*"          " "         
## 4  ( 2 )  " "          " "        "*"          " "          " "         
## 5  ( 1 )  " "          " "        "*"          "*"          " "         
## 5  ( 2 )  " "          " "        "*"          "*"          " "         
## 6  ( 1 )  " "          " "        "*"          "*"          " "         
## 6  ( 2 )  " "          " "        "*"          "*"          " "         
## 7  ( 1 )  " "          " "        "*"          "*"          "*"         
## 7  ( 2 )  " "          " "        "*"          "*"          " "         
## 8  ( 1 )  " "          " "        "*"          "*"          " "         
## 8  ( 2 )  " "          " "        "*"          "*"          "*"         
## 9  ( 1 )  " "          " "        "*"          "*"          "*"         
## 9  ( 2 )  " "          " "        "*"          "*"          "*"         
## 10  ( 1 ) " "          " "        "*"          "*"          " "         
## 10  ( 2 ) " "          " "        "*"          "*"          "*"         
## 11  ( 1 ) "*"          "*"        "*"          "*"          " "         
## 11  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 12  ( 1 ) "*"          "*"        "*"          "*"          " "         
## 12  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 13  ( 1 ) "*"          "*"        "*"          "*"          " "         
## 13  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 14  ( 1 ) "*"          "*"        "*"          "*"          "*"         
## 14  ( 2 ) "*"          "*"        "*"          "*"          "*"         
## 15  ( 1 ) "*"          "*"        "*"          "*"          " "         
## 15  ( 2 ) "*"          "*"        "*"          "*"          " "         
##           EconomicFreedom RegulatoryPolicy PersonalFreedom OverallFreedom
## 1  ( 1 )  " "             " "              " "             " "           
## 1  ( 2 )  " "             " "              " "             " "           
## 2  ( 1 )  " "             " "              " "             " "           
## 2  ( 2 )  " "             " "              " "             " "           
## 3  ( 1 )  " "             "*"              " "             " "           
## 3  ( 2 )  " "             "*"              " "             " "           
## 4  ( 1 )  " "             "*"              " "             " "           
## 4  ( 2 )  " "             "*"              " "             " "           
## 5  ( 1 )  " "             "*"              " "             " "           
## 5  ( 2 )  " "             "*"              " "             " "           
## 6  ( 1 )  " "             "*"              " "             " "           
## 6  ( 2 )  " "             "*"              " "             " "           
## 7  ( 1 )  " "             "*"              " "             " "           
## 7  ( 2 )  "*"             "*"              " "             " "           
## 8  ( 1 )  "*"             "*"              " "             " "           
## 8  ( 2 )  " "             "*"              " "             " "           
## 9  ( 1 )  "*"             " "              " "             " "           
## 9  ( 2 )  " "             "*"              " "             " "           
## 10  ( 1 ) "*"             "*"              " "             " "           
## 10  ( 2 ) "*"             " "              " "             " "           
## 11  ( 1 ) "*"             "*"              " "             " "           
## 11  ( 2 ) " "             "*"              " "             " "           
## 12  ( 1 ) "*"             "*"              " "             " "           
## 12  ( 2 ) " "             "*"              " "             " "           
## 13  ( 1 ) "*"             "*"              " "             " "           
## 13  ( 2 ) " "             "*"              " "             " "           
## 14  ( 1 ) " "             "*"              " "             "*"           
## 14  ( 2 ) " "             " "              "*"             "*"           
## 15  ( 1 ) "*"             "*"              "*"             " "           
## 15  ( 2 ) "*"             "*"              " "             "*"